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Abstract 

The first results of a new three-dimensional, finite temperature Skyrme-Hartree-Fock+BCS study 
of the properties of inhomogeneous nuclear matter at densities and temperatures leading to the 
transition to uniform nuclear matter are presented. Calculations arc carried out in a cubic box 
representing a unit cell of the locally periodic structure of the matter. A constraint is placed on 
the two independent components of the quadrupole moment of the neutron density in order to 
investigate the dependence of the total energy-density of matter on the geometry of the nuclear 
structure in the unit cell. 

This approach allows self-consistent modeling of effects such as (i) neutron drip, resulting in a 
neutron gas external to the nuclear structure, (ii) shell effects of bound and unbound nucleons, (iii) 
the variety of exotic nuclear shapes that emerge, collectively termed 'nuclear pasta' and (iv) the 
dissolution of these structures into uniform nuclear matter as density and/or temperature increase. 
In part I of this work the calculation of the properties of inhomogeneous nuclear matter in the 
core collapse of massive stars is reported. Emphasis is on exploring the effects of the numerical 
method on the results obtained; notably, the influence of the finite cell size on the nuclear shapes 
and energy-density obtained. Results for nuclear matter in beta-equilibrium in cold neutrons stars 
are subject of part II. The calculation of the band structure of unbound neutrons in neutron star 
matter, yielding thermal conductivity, specific heat and cntrainmcnt parameters, will be outlined 
in part III. 

Calculations are performed at baryon number densities of nb = 0.04 - 0.12 fm~^, a proton fraction 
of yp = 0.3 and temperatures in the range - 7.5 MeV. A wide variety of nuclear shapes are shown 
to emerge. It is suggested that thermodynamical properties change smoothly in the pasta regime 
up to the transition to uniform matter; at that transition, thermodynamic properties of the matter 
vary discontinuously, indicating a phase transition of first or second order. The calculations are 
carried out using the SkM* Skyrme parameterization; a comparison with calculations using Sly4 
at rzb = 0.08 fm^^, T = MeV is made. 

PACS numbers: 21.60.Jz,97.60.Bw,21.65.-f,21.65.Mn,26.60.Gj,26.60.Kp,26.60.-c 
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I. INTRODUCTION 



Observational properties of supernovae (SNe) and neutron stars (NSs) can serve as pow- 
erful constraints on the properties of bulk nuclear matter in density and temperature regimes 
inaccessible to laboratory experiments. In order that the theoretical description of nuclear 
matter be constrained by those properties, that description should be, as far possible, con- 
sistent over the whole extent of the relevant parameter space, a formidable task given that 
such matter spans fifteen orders of magnitude in density and reaches temperatures of up to 
lO^^K in astrophysically relevant situations. 

The main theoretical contribution to the description of nuclear matter comes in the form 
of the equation of state (EoS) of nuclear matter, obtained once the composition of the 
matter is specified. Such information is sufficient for the description of the bulk quasi- 
static properties of stars; however, in recent years increased sophistication in the modeling 
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3, |4|, of oscillations at or near 
, y, lid] have demanded a wider 



of collapsing stellar cores of the cooling of NSs 
their surface [sl, gI, and of glitches in pulsar timing 
range of physical parameters to be theoretically determined for use in hydrodynamic and 
thermal transport calculations. Such nuclear matter properties as specific heat capacity, 
electrical and thermal conductivity, shear moduli, entrainment parameters and neutrino 
opacities should be calculated within the same theoretical framework and at the same level 
of approximation as the EoS It is to this end that we begin a self-consistent study of 
the properties of nuclear matter for use in stellar modeling. 

The wide range of parameter space realized in SNe encompasses several possible phases of 
matter. Above nuclear saturation density rig ^ 2.5 x 10^^ g cm~^ nuclear matter is uniform; 
at a density a little below nuclear saturation density nuclear matter becomes unstable to 
clustering 12j and becomes inhomogeneous. The exact transition density has been shown 
to be sensitive to details of the nuclear matter EoS; at zero temperature, it is below ~ 



0.7ns 
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151 ] ■ and decreases with increasing temperature. Above a certain temperature 



T ~ 10 — 20 MeV, the inhomogeneous phase does not exist at all. At T = MeV and 
below ~ 0.2ns, the inhomogeneities are manifest as roughly spherical quasi-nuclei coexisting 
with a uniform fiuid of neutrons, electrons (ensuring charge neutrality), and a degenerate 
gas of trapped neutrinos. Equilibrium with respect to weak interaction processes is not 



reached on the timescale of collapse and the proton fraction is roughly constant ~ 0.3 16 
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17l |. As the density of nuclear matter increases, the separation of the quasi- nuclei becomes 



comparable with their physical extent: the scale on which the nuclear Coulomb interaction 
and surface tension act become similar. The matter becomes frustrated, characterized by 
a large number of local energy minima rather than a single distinct ground state. The 
transition to uniform matter through the density regime of frustrated matter is mediated by a 
series of more gradual transitions in the shape of the quasi-nuclei into more exotic structures: 



rods, slabs, cylindrical holes and bubbles, collectively termed 'nuclear pasta' Q, Q, l2Q|. 
These structures have been determined on geometrical grounds by the competition between 
the nuclear surface tension and Coulomb repulsion of adjacent nuclear formations, and 
appear to be generic properties of frustrated matter; indeed one finds interesting parallels 
with structures observed in terrestrial complex fluids 2l\ . 

The pasta phases can contribute 10-20% of the mass in the later stages of a collapsing 
stellar core [22]. As such, it is important to develop realistic, self-consistent models for this 
region of transition from inhomogeneous to uniform matter, and to determine the pasta phase 
diagram and its dependence on nuclear matter properties (symmetry energy, compressibility) 
which are somewhat uncertain at sub-saturation densities and finite temperatures 3]. 

In most theoretical models of inhomogeneous nuclear matter it is assumed that at a 
given temperature and density the matter is arranged in a periodic structure throughout 
a sufficiently large spatial region for a unit cell to be identified. Given this, the properties 
of only one unit cell need be calculated. The basic properties of the pasta phase were 



first computed using semi-classical liquid drop models or Thomas- Fermi methods [15|, |24 |. 
The energetically preferred pasta phase is obtained by calculating the energy-densities of the 
separate phases (including that of uniform matter) across the whole density and temperature 
regime in which pasta is expected to exist, and selecting the phase that gives the lowest 
energy-density at a given baryon number density. This means that the nuclear shapes 
expected to appear must be specified a priori. 

Fully microscopic one-dimensional (ID) Hartree-Fock (HF) calculations of SN matter 
were carried out by Bonche and Vautherin 25|, l26|] (from hereon B V) . These calculations did 
not include a systematic study of possible spurious effects, e.g. shell effects caused by the 
discretization of the unbound neutron energy spectrum by the finite computational domain. 

While the ID-HF method self-consistently describes the nuclear bulk, surface and shell 
effects, it restricts the nuclear shapes to being spherically symmetric. It requires the spherical 
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Wigner-Seitz (WS) approximation, in which the generally non-spherical unit cell is replaced 
by a spherical one with the same volume. The WS approximation is good as long as the 
nuclear structures that form the lattice points are sufficiently widely spaced. This condition 
is satisfied everywhere except close to the transition density, when their spacings become 
comparable with their individual physical extent. Particularly, the transition to uniform 
matter cannot be treated in the WS approximation; the spherical unit cells by their nature 
do not continuously fill the space. In order to treat the transition to uniform matter correctly, 
three-dimensional (3D) studies are necessary. 



The techniques of quantum molecular dynamics (QMD) and similar [27|] are semi-classical 
microscopic approaches that lead to an improved treatment of the pasta regime. In this 
method a large number of nucleons is dynamically evolved in a large cubic box with periodic 
boundary conditions and without assumption on the nuclear shape. The box can be made 
large enough to include the effects of electron screening on the Coulomb lattice energy. 
The pasta's dynamical response to a neutrino flux has been investigated in this framework 
28l . l29l | , and signiflcant strength at low energies from excitations of the internal degrees of 
freedom of the pasta has been found. This may be an important effect in the treatment 
of neutrino interactions with matter in SN simulations, especially concerning shock revival. 
The pasta shapes themselves and their sequence have also been studied in detail using QMD 



at zero and flnite temperature 
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31| . It has been shown that more complicated shapes 



intermediate between the canonical pasta shapes may exist, and that pasta can also be 
formed by compression of a bcc lattice of spherical nuclei on a timescale much shorter than 
that of core collapse. The drawback of the QMD method is that the effective interaction 
used is very schematic and that the important microscopic shell effects are not included. 

Recently, the three-dimensional (3D) HF method in a cubic box using the Sly4 Skyrme 
parameterization was applied to neutron star matter by Magierski and Heenen, and Gogelein 



and Muther 
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331]. Pasta shapes coexisting with the external neutron gas emerged natu- 



rally, and in addition other exotic shapes were found because no prior selection of the nuclear 
shapes expected was required. Magierski and Heenen examined the shell structure of the 
external neutron gas, which, it was posited, arises from the scattering of the free neutrons 
comprising the external neutron gas off the nuclear clusters. This effect has two related 
consequences: (i) the energy distribution of the free neutrons is discretized like that of the 
bound nucleons, forming a shell structure; (ii) the scattering causes an effective interaction 
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between nuclear clusters analogous to the Casimir effect in quantum field theory (and led it 
to be dubbed the Fermionic Casimir effect) 3J]. The energy associated with that interac- 
tion is comparable to the energy difference between the different shape phases. As a result, 
the order in which the nuclear shape changes occur with changing density may be different 
to the sequence postulated in simpler models, and several shapes may coexist at the same 
density in different areas of the star, destroying the long range periodicity of the matter. 

Both 3D-HF studies to date have calculated nuclear configurations at a limited number 
of values of densities and number of nucleons in the unit cell, and only for proton fractions 
expected to be found in neutron star crustal matter. Only Gogelein and Muther [s^ per- 
formed calculations at finite temperature. To self-consistently probe the energy of various 
pasta shapes, both independent quadrupole moments (g2o and ^22) should be constrained; 



the study of Magierski and Heenen 



321 ] imposed a constraint only on the ^20 component 



of the proton quadrupole moment, whereas Gogelein and Muther performed unconstrained 
calculations. Finally, it is not clear what role the numerical method itself plays in the results 
obtained; for example, what are the effects of the finite computational cell? 

In the present work the first comprehensive 3D-HF study at finite temperature with a 



Skyrme-type interaction 
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38| of the properties of inhomogeneous nuclear matter 



near the transition density to uniform matter in SN matter is presented. Preliminary results 
were outlined earlier jsol; in this paper a full account of the model is given including an 
extensive survey of the computational details necessary to prove validity of the approach 
to modeling nuclear pasta. Special attention is paid to the role played by the finite size of 
the computational cell. Furthermore, this work concentrates on a systematic exploration of 
the 'shape phase space' of nuclear pasta using a strict constraint on the two independent 
quadrupole moments of the neutron density distribution to probe the energy of various 
nuclear configurations. The 3D-HF calculation across the phase transition region allows the 
EoS of nuclear matter to be developed self-consistently from neutron drip density (where a 
ID-HF model is sufficient) right through to uniform matter densities, avoiding any artificial 
discontinuities that might occur matching together EoSs covering different density regions 
which has been a common practice so far. 

The computational scheme presented in this paper is constructed to be flexible enough 
to easily adopt a variety of nuclear interactions other than those of the Skyrme type, and 
to study different boundary conditions. Its application to beta-equilibrium neutron star 
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matter will be presented in part II of this work. Inclusion of more general Bloch boundary 
conditions and results of the first self-consistent calculation of the band structure of neutrons 
in the pasta phases of neutron star matter will be presented in the third paper (part III). 

The paper is organized as follows: the theoretical method is described briefly in Sec. [ITl 
Implementation of the model in the form of the TAMAR code is outlined in Sec. Ill A[ Nu- 
merical application to terrestrial nuclei and to inhomogeneous matter is detailed in Sees. Ill Bl 
and III CI In Sec. Ulllthe results for supernova matter are presented and discussed. Conclu- 
sions and outline of future developments are given in Sec. [IVl 



II. THE SKYRME-HARTREE-FOCK (SHE) + BCS MODEL 



The computational framework used is briefly outlined in this section. Details relevant to 
the calculation of inhomogeneous nuclear matter in stars are given; more complete descrip- 
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tions of the Hartree-Fock method are available elsewhere 

In the HF method, the ground state wavefunction is approximated by a single Slater 
determinant |$) = \ni,n2, ) = a\a2 |0) which is found by minimization of the expec- 
tation value of the hamiltonian of the system S{^\H\^) = 5£^skyrme ['^'] = 0, where £^skyrme is 
the energy- density functional 



381 ]. Minimization with respect to the single particle wave- 
functions \ni) yields a set of single-particle Schrodinger equation with the one-body HF 
potentials Ug-. 



^ ^ /X / X (V X a) 



(1) 



Here, q = p,n labels the isospin states, i the single particle states, U(so),g is the spin-orbit 
potential, and m* is the effective mass. 



Using the Skyrme interaction the one body potentials are obtained 
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37j: 



1 



(2 + «)(! + ix3)p - 2(1 + X3)pg - ail + 



pI + pI 

P 



- it4(V ■ J + V ■ Jg) , 



(2) 
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U(so),g = ^^4(Vp + Vpq) + ^(tl - t2)Jq " ^(a^l^l + a;2t2)J, (3) 

and the effective mass is 



+ -[ti(l + + t2(l + \X2)\P - -[tl{\ + Xi) - t2{\ + X2)]pg . (4) 



where and Xi are parameters of the Skyrme interaction, p = Pp + Pn are the nucleon 
densities, r = Tp + r^ are the kinetic energy densities and J = Jp + Jn is the spin current. Two 



sets of Skyrme parameters were used in the present work; the SkM* parameterization [43 1 
based on a high precision description of nuclear ground states as well as surface energies and 
fission barriers of nuclei is used throughout the paper, and the Sly4 parameterization 4^ 
developed to describe neutron rich nuclei and neutron matter and biased more towards 
astrophysical applications. 

In order to reduce the computational task in the present work the spin-orbit force has 
been not included in the HF Hamiltonian. Undoubtedly the spin-orbit interaction is an 
important component in determining the correct shell energy and single particle spectrum. 
However, the initial aim is to set up and thoroughly test the computational procedure, and 
to examine the optimum way to apply it to the calculation of the sub-nuclear EoS. The 
results are to be seen as a reference set with which later developments can be compared, 
and, in this case, the effect of the spin-orbit force systematically examined. 

BCS pairing correlations are added to the Skyrme energy-density functional at zero tem- 
perature as outlined in [42]. At finite temperature the single particle level occupation prob- 
abilities are given by the Fermi-Dirac distribution. 

In the uniform matter case single particle wavefunctions can be represented by plane 
waves, the gradient terms in Eq. ([2]) disappear and the expressions for the single particle 
potentials and the total energy- density of nuclear matter are analytic. A full explanation 
of obtaining the uniform matter EoS using the Skyrme-Hartree-Fock model can be found 
in 4^. 



A. Computational Implementation 

Eqs. ([1]) are solved in co-ordinate space in one cubic unit cell of the locally periodic nuclear 
matter. Note that this choice admits several different crystal lattice types; simple cubic (sc), 
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body centered cubic (bcc) and face centered cubic (fee). To reduce the computational task 
further, only nuclear configurations that conserve reflection symmetry in the three Cartesian 
directions are considered. It follows that the computation need be performed only in one 
octant of the unit cell. 

It is expected that the absolute minimum energy of the cell is not going to be particularly 
pronounced and there will be a host of local minima separated by relatively small energy 
differences corresponding to different nuclear geometries. In order to systematically survey 
the 'shape space' of all nuclear configurations of interest, the deformation of the nuclear 
configuration must be controlled via a constraint imposed on the nuclear density distribu- 
tion. Reflection symmetry across the three Cartesian axes is assumed in the present model, 
thus eliminating asymmetric deformations such as dipole or octupole. A constraint is thus 
imposed on the quadrupole moment of the nuclear densities; the next order deformation 
consistent with the shape of the unit cell used, hexadecapole, is expected to give energy 
variations an order of magnitude smaller than those of the quadrupole deformation. The 
constraint is applied only to neutrons; because the matter that is modeled is neutron rich, 
the proton distribution is expected to follow the neutron distribution to a good degree of 
accuracy. The scheme used to implement the constraint can be found in |461] . Specific defor- 
mations will be referred to in terms of the standard 'polar' deformation co-ordinates 7) 



[see, e.g. 



47|). 



The computational grid has a spacing Ax, Ay and A;z and is arranged so that the ith 
collocation point is located at Xi = {i + |) Ax and similarly for the y and z directions, where 
i is an integer. The boundary conditions are simple periodic; a given function in the cell 
must obey 0(r + T) = 0(r). where T is the translation vector from the position r to 
the equivalent positions in the adjacent cells. These conditions are enforced by calculating 
derivatives and solving the Poisson equation for the Coulomb potential on the grid 4^ in 
Fourier space. The Coulomb solver utilises the FFTW software package In NS matter. 



the more general Bloch boundary conditions must be used 49|, |50|, but in SN matter (finite 



temperature and small density of dripped neutrons) this is not necessary. The Coulomb 



51| . Integrals of functions are 



exchange term is evaluated via the Slater approximation 
calculated using the trapezoidal rule; for periodic functions this is exact [52l |. 
Two iteration methods were used in this work. The imaginary time step 



531 was used 



for the first tens of iterations for its robustness, and the damped gradient step |46|, |5J] used 



9 



afterwards for its quicker convergence at later iterations. The iteration is considered to 
have converged to the desired accuracy once the total variance in proton and neutron single 
particle energies (A/i^)^ = Y^i^hg (0i,<?l(^HF)^l0i,g) " ((^i.gl^Hpl^i,?))^ becomes less than 
1.0 keV^ which ensures convergence in the total energy to less than one part in 10®. 

A uniform background of electrons is included to ensure charge neutrality. The electrons 
generate their own Coulomb potential, $e which is determined in the same way as the proton 
Coulomb potential. Screening of the electron Coulomb potential is taken into account simply 
by omitting the infinite part of the potential. Free Fermi gas expressions for the energy and 
entropy of the electrons are used 42|. The electrostatic interaction between the electrons 
and the inhomogeneous proton distribution is taken into account by adding to the single 



25| 



particle kinetic energies of the electrons an electrostatic potential term as in 

The total energy of the nuclear configuration including the contribution of the Coulomb 
force from the electrons and protons is given by 37l |: 



^skyrme ~t" ^coul 2 ( ^kin ~l~ ^ ^ p(^fS j "l~ i^roarr ~\~ ^^quad ~l~ ^^pair ~l~ 
^ /3 ^ 



exch,coul 



~\~ ^e.coul ~l~ ^1 



lattice ) 



(5) 



where the various contributions are, respectively, from the kinetic energy of the nucleons, the 
single particle energies, the Skyrme rearrangement energy, the energy due to the quadrupole 
constraint, the pairing energy, the Coulomb exchange energy, the Coulomb interaction of 
the electrons and the lattice energy. 

At finite temperature the Helmholtz free energy F = Stot — TS is calculated. Here T 
is the temperature and S is the total entropy of all particles present calculated using the 
free Fermi gas expressions derived in Appendix A of Ref. 4^]. The free energy- density is 
then obtained by dividing by the volume of the unit cell / = F/V. The pressure density of 
the matter is calculated using p = — / + ^ . /ijnj where the sum is over all particle species 
(labeled i), fii are the chemical potentials and the number densities. The evaluation of 
the chemical potentials is given in detail in 42 1. 

As a final note, in the following the mean number densities of particles in the matter 
(that is, their bulk values) are simply denoted nn, n^, etc. Where necessary, the local number 
densities will be given as a function of position in the unit cell nn{x, y, z), etc. 
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B. Calculation of terrestrial nuclei 



The calculation scheme described above has formed the basis of the TAMAR computer 



code. The code has been extensively tested 



42j for terrestrial nuclei against two reference 



HF codes. The first is a 2D code SKYAX in cylindrical polar co-ordinates which enforces 
axial symmetry 55|]. This code uses a finite difference scheme to treat the derivatives 
and calculates the Coulomb potentials using a separate iterative process. The second is a 
three-dimensional code TDHF3D [sgI. Like the TAMAR code, it represents derivatives and 
solves the Coulomb force using Fourier Transforms, allowing for a more direct comparison. 
Both reference codes use the damped gradient iteration step. Since the three codes use 
different prescriptions for the calculation of the Coulomb and pairing forces, the test has 
been performed with the basic Skyrme nuclear interaction without the constraints or the 
Coulomb and pairing forces. Two Skyrme parameterizations SkM* and Sly4 have been 
used to calculate the binding energies and single particle energies of a selection of doubly 
magic nuclei and ^^Fe and excellent agreement with both reference codes was found. As an 
example, in Table [T] the binding energies and root mean square proton and neutron radii as 
calculated by the three codes using the Sly4 Skyrme parameterization are compared. 

Furthermore, the method used in the TAMAR code to calculate the Coulomb potentials 
has been carefully examined. This method is intrinsically periodic; when solving for a nucleus 
the presence of a surrounding simple cubic lattice of the same nuclei is implicitly taken into 
account. There is thus a contribution to the binding energy of the nucleus from the lattice. 
The lattice energy decreases with increasing the computational volume; for large volumes, 
effects of the finite size of the nucleus can be neglected and the lattice energy is given by 
that of a simple cubic lattice of point like nuclei of charge Ze\ El = —1.444-^. Testing of 
the variation in total energy with cell size revealed that the Coulomb energy tends towards 
a constant value at large a as expected, and the contribution of the lattice energy falls off 
as The proton Coulomb energy was found to be in agreement with the TDHF3D 

reference code. 
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C. Inhomogeneous matter 



1. General Numerical Considerations 

For each run of the code at a given temperature T and baryon number density nt = A/V, 
where A = N + Z is the number of nucleons (A^ neutrons and Z protons), and V the volume 
of the unit cell, up to four free parameters need to be specified: A, the proton fraction = 
Z/A, and the quadrupole moment parameters (3 and 7. Clearly at a given value of rib, A and 
V are not uniquely determined and calculations have to be performed over a range of values 
of A, or equivalently V, at that given density to find the value which gives the minimum 
free e„ergy-de,.,t. The p.oto. fraction .s kept constant tlnoughout th,s work at ^ 0^ 
a value appropriate for collapsing stellar cores [15|. If the quadrupole constraint is applied, 
the deformation parameters j3 and 7 are also free parameters and must be varied over the 
deformation space, searching for shapes corresponding to minima of the total energy. 

To build an EoS of the pasta phases of inhomogeneous matter, such a series of calcula- 
tions has to be performed at each point of the parameter space (nb, T, i/p) with adequate 
spacing, which is extremely computationally intensive. Part of the task will be to find the 
optimum strategy to obtain adequate coverage of parameter space; for example, some areas 
of parameter space are physically less important than others; this will lead to selection of 
those values of the parameters that are most likely to be physically manifest, and a corre- 
sponding reduction in the computation time. Examples of this selection process are given 
in Sec. HH 

The calculation of inhomogeneous nuclear matter requires consideration of some numeri- 
cal effects that are common to the calculation of isolated terrestrial nuclei and, in addition, 
a number of numerical effects specific to the particular problem at hand. The computational 
volume V is uniquely defined by the baryon number density nb and the number of nucleons 
A present in that volume by = A/n\^. Since V is fixed at a given nb, V^, the number 
of grid points and the grid spacing are inversely related - a constraint not present in the 
calculation of isolated nuclei. There is a balance between the number of grid points used 
to obtain the desired accuracy and the required computational time. The number of grid 
points determines the maximum number of orthogonal wavefunctions that can be repre- 
sented on the grid. For the grid parameter of interest in this work {Axa ~ 1 fm, rix^ ~ 10), 
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this number is in the thousands. Ideally the calculation should be performed using all the 
available wavef unctions as it is not known a priori which is the preferred set of occupied 
states forming the configuration corresponding to the minimum total energy, and yielding a 
realistic approximation to the physical ground state. Furthermore, with increasing tempera- 
ture a larger number of single particle states will become occupied, so a smaller value of grid 
spacing will be required for higher temperatures. However, in the interests of reducing the 
computation time, the number of wavefunctions participating in the iteration process should 
be minimal. As a compromise, each calculation begins by including all the wavefunctions 
available for a given grid; then, after a certain number of iterations, those wavefunctions 
which have occupation amplitude below 10^® are excluded from the iteration process after 
each iteration. It was found that 100 iterations was the safe choice of the cut-off point. 

The optimal grid spacing was determined as follows. Calculations were performed at 
a proton fraction of = 0.3, baryon number densities of nb = 0.06, 0.08 and 0.10 fm~^, 
temperatures T=0, 2.5, 5.0 and 7.5 MeV, quadrupole deformations 7 = 0° and P = 0.0 and 
1.0, and nucleon number A = 100-1200 with a step of 20. For each of these collections of 
points in parameter space, an average grid spacing of (Axa) = 1.0, 1.1, 1.2 and 1.3 fm was 



tested 42j, and the difference in energy obtained at each grid spacing was calculated. It 
was found that using a grid spacing of 1.0 fm guarantees an accuracy 1 part in 10"'^ in the 
energy for temperatures up to 5 MeV. At 7.5 MeV, a smaller grid spacing is required. The 
calculations presented in this paper do not use grid spacings below 1.0 fm at this temperature 
because of the large computation time required, so care has to be taken in interpreting the 
results at this high a temperature, where, in principle, higher momentum states need to be 
included. 

Another important requirement the model has to satisfy is the independence of the final 
outcome on the form of the initial wavefunctions as mentioned in Sec. ??. Two cases were 
considered: (i) both proton and neutron wavefunctions are of the form Gaussian times 
polynomial (to assure orthogonality), (ii) the proton wavefunctions are Gaussian and the 
neutron wavefunctions are plane waves. It was found that both cases converged to the same 
energy to less than 1 part in 10~^. 
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2. Effect of Finite Cell Size on the Computation of Inhomogeneous Matter 

When performing calculations of isolated nuclei the computational volume can be made 
much larger than the size of the nucleus, and is constrained only by the computation time. 
It follows that the effect of the choice of the computational volume can be minimized. In 
contrast, the nuclear configurations in stars span the entire computational volume; indeed 
the volume itself acquires a physical meaning. Thus the numerical effects of the size of a 
finite cell must be carefully examined. 

Given a particular nuclear shape obtained in a cubic cell, it should in principle be possible 
to double the cell size and, searching over deformation space, find a configuration equivalent 
to two instances of the original shape adjacent to each other. The free energy- density 
obtained in the doubled cell should be identical to that in a single cell. 

To examine this effect, calculations were performed in a box that is double the length 
in the z-direction compared to the other two directions. The expectation was to find a 
configuration similar to the one in a cubic box at the same baryon number density (with 
half the number of nucleons), with the same free energy energy. Fig. [T] displays one of such 
configurations. The nucleon density distribution has been rendered in three dimensions, 
with blue indicating the lowest densities and red the highest. The right-hand cell displays 
the nuclear shape obtained within a cubic cell; it is similar to the 'lasagna' phase, with an 
additional cylindrical bridge joining the slabs. The left-hand cell displays a shape that is 
very similar to two instances of the modified lasagna shape adjacent to each other. The free 
energy densities of the two configurations agree to within 1 part in 10*^. 

Fig. [T] demonstrates clearly that the nuclear structures obtained are not artifacts of the 
finite cell size (and more evidence will be presented in paper II for neutron star matter). 
However, the finite cell size can be a source of spurious shell effects. In analogy with a Fermi 
gas in a box, shell effects caused by the discretization of the physical space due to the finite 
computational volume may occur. These effects will manifest themselves especially at high 
densities and temperatures when a large number the nucleons are unbound. To illustrate the 
form of these numerical shell effects, the free energy-density is plotted in Fig. [2] as a function 
of A in a wide region up to A=2000 at a density of ny, = 0.11 fm~'^, (/?, 7) = (0.0,0°) and T = 
2.5 MeV. The pronounced oscillations in the free energy-density curve are entirely numerical 
in origin. With increasing box size and temperature the bulk energy- density of the matter 
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in the box approaches that derived from the semi- analytic formulae for the free energy- 
density of a free nucleon gas. These numerically induced shell effects manifest themselves 
in a form distinct from the physical shell effects, arising due to a combination of the shell 
energies of bound nucleons and unbound neutrons scattered by the bound nucleons, which 
are characterized by more rapid fluctuations in nucleon number A typically at lower densities 
and temperatures and lower values of A. The distinction between the form and occurrence 
of the two types of shell effects is encouraging as it allows their easy identification. In 
such a situation where the shell effects are purely spurious, the physical value of the free 
energy-density is not the minimum, but that value to which the free energy tends at high 
A. 

In order to accurately locate free energy-density minima as a function of A and in energy- 
density surfaces, the numerical shell effects, which introduce spurious minima, should be 
identified and corrected for. 

III. RESULTS AND DISCUSSION 

Keeping the proton fraction fixed at yp = 0.3, and at selected densities and temperatures 
nb and T, the variation of free energy- density of inhomogeneous nuclear matter / with the 
parameters A, (3 and 7 is now explored. 

A. Variation of free energy-density with nucleon number at constant deformation 

The minimization of free energy-density with respect to total nucleon number A (equiv- 
alent to unit cell size) at constant deformation is illustrated in Fig. [3] at two densities rib = 
0.06 fm~^ (top two plots) and 0.10 fm~'^ (bottom two plots). The deformation is held con- 
stant at (/5, 7) = (0.0,0°) for the left two plots and (/5,7) = (1.0,0°) for the right two plots. 
Each plot displays results obtained T = MeV and T = 2.5 MeV to demonstrate the effect 
of temperature on the form of the curves. Note that in general, the lower the temperature, 
the higher the free energy- density: increasing the temperature means the amount of useful 
energy (i.e. available for mechanical work) decreases as the entropy increases. In addition, 
the free energy-density versus A at rib = 0.06 fm~^, T = 5 MeV, and (/3, 7) = (0.0,0°) is 
shown in Fig. HI The variation of important EoS quantities (free energy- density /, entropy 
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density s and pressure density p) with A is tabulated in Table HTl for nb = 0.06 fm ^, T = 
2.5 MeV and (/5,7) = (1.0,0°). 

The form of the free energy curves displayed is effectively a superposition of physical 
effects related to the geometry of the nucleon density distribution, physical shell effects and 
numerical shell effects. There are also occasional discontinuities in the curves (e.g. at A = 
820 in the plot at = 0.06 fm'^, T = MeV and (/3,7) = (0.0,0°)) caused by a sudden 
change in the geometrical configuration of the nucleons in the unit cell, as will be discussed 
below. 

The competition between nuclear surface tension and Coulomb energy manifests itself as 
a broad minimum. Because such minima have their origin in the nuclear geometry, they 
will be referred to from now on as geometrical minima. It is upon this broad curve that 
the shell effects are superimposed. The physical shell effects are particularly evident at 
zero temperature as irregular oscillations in energy which are responsible for local minima 
appearing in the curves. As noted in 25|], in matter with a proton fraction of 0.3, shell 
effects yield minima far from the familiar magic numbers of terrestrial nuclei. Shell effects 
have their highest magnitude at low values of A. This is because the density of states is 
lowest at smallest A; the energy gaps between shells are at their largest, thus the difference 
in total energy between a filled shell and, say, a shell occupied by only two nucleons, is 
relatively large. With increasing A, the density of states increases, the shell separations 
become smaller, and the effects on total energy diminish. 

As the temperature increases, the shell effects become less pronounced, disappearing 
when the typical spacing between shells is less than k^T. This is particularly evident as 
one follows the form of the free energy curve at nb = 0.06 fm~'^, = (0.0,0°) as one 

increases temperature from T = MeV through T = 2.5 MeV to T = 5 MeV (top left plot 
of Fig. [3] and Fig. Hj). At T = 5 MeV the shell effects have completely disappeared and 
only the broad geometrical minimum remains aX A = 600. The plots of free energy-density 
with increasing A at rib = 0.06 fm~^ and zero deformation is consistent that found in the 
ID-HF calculations of Bonche and Vautherin 25|]. The absence of the spin-orbit interaction 
alters the details of shell effects but does not change the picture qualitatively. The broad 
minimum appears at A = 600 in both works. 

As deformation increases to (/3,7) = (1.0,0°) at the lower density rib = 0.06 fm~^, the 
curves remain qualitatively similar. The shell effects at T = MeV have increased in 
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magnitude a little; indeed, the minimum appears at A = 300 due to shell effects rather 
than the nuclear geometry. The geometrical minimum occurs at around the same value of A 
despite the different nuclear geometry, although this is not always the case. The geometry at 
7) = (1.0, 0°) is displayed in Fig. HJ and the geometry at 7) = (0.0, 0°) in the top left 
panel of Fig. [T^l The free energy- density at the geometrical minimum is about 2 keV fm^^ 
greater at (/?,7) = (1.0,0°) than at zero deformation. Several EoS quantities are tabulated 
in Table [IT] as they vary with A at rib = 0.06 fm'^, T = 2.5 MeV, (/5, 7) = (1.0, 0°). One can 
see that although the free energy changes relatively smoothly with A at this temperature 
as most of the shell effects have been washed out, quantities such as pressure and entropy 
which depend on the derivative of the free energy still fluctuate a significant amount as A 
increases. 

The bottom two plots of Fig. [3] show the variation of free energy-density with A at a 
higher density nt, = 0.10 fm"^. At zero deformation the form of the curves is qualitatively 
different from those at the lower density. The nucleon distribution is close to being uniform 
at this density (in fact there are small bubbles present in otherwise uniform nuclear matter; 
see the top two panels of Fig. [T3l) and the shell effects are dominated by those of a uniform 
Fermi gas in a cubic box as discussed in Sec. Ill C 21 and illustrated for completely uniform 
matter in Fig. [2l As has been discussed, these shell effects are numerical in origin and 
show a characteristic form significantly different to physical shell effects in bound nuclear 
structures; they are more regular and persist to higher values of A and T. At this density, 
since structure is still present, they are overlaid on a broad geometrical minimum. The 
spurious shell effects give local minima deeper than the geometrical minima, so until they 
are subtracted off, it is no longer good practise to take the ground state of matter to be 
that corresponding to the minimum energy configuration. At zero temperature there are 
still some physical shell effects visible; these become washed out at T = 2.5 MeV leaving 
behind the purely spurious effects superimposed on the geometrical minimum. 

At rib = 0.10 fm~'^, the free energy-density of significantly deformed matter {j3 = 1.0) is 
appreciably higher than that of non-deformed matter. This is an indication that the transi- 
tion to uniform matter is being approached: more uniform configurations are energetically 
preferred rather than highly deformed ones. The variation of / with A has a form similar 
to that at rib = 0.10 fm~^; the matter, being deformed, no longer exhibits the spurious shell 
effects associated with uniform nuclear matter. 
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3D renderings of a selection of neutron density distributions taken at T = 0.0 MeV and 
T = 5.0 MeV for the sequences of free energy- density versus A at = 0.06 fm~^, (/3,7) 
= (0.0, 0°) (see top left panel of Fig. [3] and Fig. H]) are shown in Fig. O The unit cells are 
displayed to scale relative to each other, so the increase in cell size as A increases is clearly 
visible. The smallest cell width is 11.9 fm and the largest is 27.1 fm. It can be seen that 
the constraint on both quadrupole moments has ensured that, for the most part, the same 
geometry is obtained for each value of A, at constant {[3, 7). Thus a unique geometrical 
configuration could be followed as the total nucleon number A is increasing. The geometry 
also remains the same for the two temperatures shown. At zero temperature there is no 
neutron gas external to the nuclear structure, although a neutron halo is visible; at T = 5.0 
MeV an external neutron gas exists. The main structure obtained is notably not a member 
of the canonical pasta collection: its shape is similar to that of three cylinders intersecting 
each other at right angles. The proton distribution, not displayed, follows the neutron 
distribution very closely. The external neutron gas is very dilute at both temperatures, as 
has been found in previous studies of SN matter 25|. 

Although for the most part, the same geometrical configuration is obtained at a given 
(/3, 7), there are some configurations that do not follow the pattern (see the bottom panels). 
This is a result of the freedom inherent in the TAMAR code from not constraining higher- 
order deformations; for example, a configuration may have the same quadrupole moment, 
but a different hexadecapole moment. These departures from the energetically preferred 
geometrical configuration are rare, however, and are clearly distinguished from others in the 
energy plots as discontinuities in the curves of order of 1 - 5 keV fm~^. One such point is 
observed at T = 0.0 MeV at a nucleon number of A = 820 and three at T = 5.0 MeV: A 
= 1080, 1140, 1200. The former has a small number of nucleons split off from the main 
structure, forming a 'mini-nucleus' arranged in a body-centered cubic lattice with respect 
to the larger structure. The latter two points shows a more significant change: now the 
nucleons have arranged themselves more equally in a body-centered cubic structure with 
nucleon 'arms' joining the nuclei at the lattice points. 
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B. Variation of free energy-density with nucleon number and deformation param- 
eters P, 7 



In the previous section the variation of free energy-density / with nucleon number was 
considered at two specific values of the deformation parameters. In this section the study 
of the variation of / is extended to include the whole range of values of (3 of interest and for 
selected values of 7. This will give a series of energy-deformation surfaces / = /(A,/?, 7 = 
const.) which contain information about shell and geometrical effects associated with varying 
nuclear shapes. 

Fig. [H] displays the energy-deformation surfaces obtained at the density values Ub = 0.08 
fm~^ and 7 = 60° for four values of temperature T: 0, 2.5, 5.0 and 7.5 MeV. It can be seen 
that the energy-deformation surfaces retain the principal features, discussed in the previous 
section, along the curves of constant deformation. The physical shell effects associated with 
bound nucleons are still visible at lower temperatures. Deformation of nuclear structure 
breaks some of the degeneracy of the single particle energies within particular shells, leading 
to a larger number of sub-shells; this is seen on the energy- deformation surfaces as the 
increasing prominence of shell effects as (3 is increased. The result of the shell effects is 
many local minima in the surfaces at zero or low temperature, of order of 5 keV fm~^ deep. 
The variation in shell energy with respect to P at constant A is expected to be caused by the 



varying shell structure of the bound nucleons. Magierski and Heenen [32|] suggest that in 
neutron star matter, the shell effects of the unbound neutrons (caused by their scattering off 
the nuclear structure) contribute significantly to the variation of energy with deformation; 
however, this is not expected to be the case here because the density of unbound neutrons is 
an order of magnitude smaller in SN matter than in NS matter at a corresponding density. 

The shell effects become washed out at higher temperature leaving behind the geometrical 
minima. We see that these occur not only with respect to A, but also with respect to (3, and 
result in several local geometrical minima appearing on each energy- deformation surface. 
For example, at T = 5.0 MeV, rib = 0.08 fm~^ and 7 = 60° (Fig. [6l bottom left) there 
are two broad minima that appear, centered at /5 ~ 0.3, A ^ 800 and (3 ~ 0.7, A ^ 800. 
The geometrical configurations they correspond to are displayed in Fig. [71 The former 
corresponds to slab like nuclear shape ('lasagna'); the latter to a cylindrical hole nuclear 
shapes. They are separated by a 'ridge' in the surface, running along j3 = 0.6 which delineates 
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the regions on the energy surface where the two geometries occur. The free energy-density 
changes smoothly in the vicinity of the minima, and the nuclear geometry locally remains 
the same. 

Fig. [H] shows the energy-deformation surfaces obtained at the density values rib = 0.06 
fm~^ and temperature T = 5.0 MeV for three different values of the deformation direction 7: 
0°, 30° and 60°. The minimum with respect to A occurs at 600-800 and is fairly constant with 
respect to (3. As 7 changes, so do the geometries obtained as jS increases. A deep minimum 
appears on the 7 = 60° energy surface at A ^ 600, P = 1.1, and indeed this corresponds 
to the absolute minimum found at this density and temperature. It is important to note, 
however, that the energy separations between local minima is relatively small (of order keV) 
and one must expect multiple configurations to exist at a given density and temperature. 

Fig. [9] shows the energy-deformation surfaces obtained at the densities n^ = 0.08 - 0.10 
fm~^ for 7 = 0° and temperature T = 2.5 MeV. As the density increases the nature of the 
shell effects is seen to change, becoming dominated by a spurious component the highest 
density and low deformations. The free energy-density rises increasingly sharply for (3 > OA 
as density increases. This is a clear indication that the favored ground state is becoming 
closer to uniform matter; it is costing more energy to deform the matter into structure. It 
is also interesting to note that the spurious shell effects disappear at large deformation, as 
the structure becomes far from uniform again. 

C. The EoS and the Transition to Uniform Matter 

By taking the minimum free energy-density configurations from each baryon number 
density, and temperature, the EoS of non-uniform matter can be constructed. We have 
performed this minimization over A, (3 and three values of 7 (0°, 30°, 60°) at seven different 
values of density: = 0.04, 0.06, 0.08, 0.09, 0.10, 0.11, 0.12 fm'^ and three different values 
of temperature T = 0, 2.5 and 5 MeV. One should bear in mind that spurious shell effects 
might alter the position of the minima in phase space at densities and temperatures close 
the transition to uniform matter. 

In Tables UTTl to IIVI the variation of relevant EoS quantities with density and temperature 
are given. Values of A, (3 and 7 corresponding to the minima are tabulated, together with 
the corresponding values of free energy-density /, entropy density s, pressure p, proton. 
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neutron and electron chemical potentials /ip, /in, /Xe, and the quantity S = r^Ms/^c, where 
'"RMS is the nucleon root mean square radius and rc is equal to half the width of the cell. S 
is a measure of the fraction of the cell occupied by the nuclear structure; for uniform matter, 
5 = 1. 

The energy minima appear at values of A that are generally smaller than those found by 
BV [25] (e.g. at T = 5 MeV, nb = 0.06 fm~^, a minimum at A = 600 is found as opposed to 
A = 900 for BV). This is a result of allowing deformed nuclear structure; similar values of 
A to BV are obtained if one only admits zero deformation configurations (it was observed 
in Sec. [TllA]that at rib = 0.06 fm'^, T = - 2.5 MeV, (/3,7) = (0.0,0°), the same value of 
A as BV for the minimum energy was obtained). 

In Fig. [10] the neutron density profiles are shown along a line connecting two adjacent 
unit cells in three perpendicular directions for the minimum energy configurations at all 
seven densities and T = 2.5 MeV. The axes are labeled x,y and z, but one should note 
that the choice of axis labeling is arbitrary in relation to the macrophysical structure of the 
matter. In Fig. [TT] 3D renderings of those neutron density distributions in one unit cell are 
displayed. It can be seen that the transition proceeds to uniform matter proceeds along the 
recognised ordering of the pasta shapes: cylindrical, planar, cylindrical hole and spherical 
hole. We finally achieve a self-consistent transition to uniform matter from increasing the 
density at around ^b = 0.10 - 0.11 fm~'^ at this temperature. 

In Fig. [12] 3D renderings of the neutron density distribution of energy-minimum config- 
urations taken at nb = 0.10 fm~'^, for T = 0, 2.5, 5, and 7.5 MeV are displayed. Up to T 
= 5 MeV, the configuration is bubble-like; at T = 7.5 MeV a self-consistent transition to 
uniform matter has been reached by increasing the temperature. 

Despite the fact that the canonical pasta shapes were obtained at the absolute minima 
in energy at the selected values of density, it should be noted that a far richer variety of 
nuclear shapes were obtained here at local minima close to the absolute ground state. It is 
expected that, as calculations are performed over a greater number of densities in between 
those already studied, the transition between the canonical pasta shapes will be mediated 
by a number of different nuclear geometries. In addition, it is likely that at a single given 
density and temperature, one will find many pasta shapes co-existing. The implication is 
that at a macrophysical level, the bulk properties of matter will be given by an averaging over 
the properties of those individual pasta phases, and that the thermodynamical properties of 



21 



matter will change smoothly with density and temperature through the pasta regime until 
the phase transition to uniform matter is reached. 

For the present model to produce the transition to uniform matter consistently, the var- 
ious thermodynamic quantities calculated should tend in some way towards those calculate 
in the uniform matter where the Skyrme energy-density is given analytically. In order to 
check this, the difference in free energy-density, entropy density and pressure between the 
inhomogeneous matter and the uniform matter Skyrme-HF model is plotted in Figs. [13] 
to [15] for the densities already considered and temperatures of MeV, 2.5 MeV and 5 MeV. 
The difference in free energy-density A///uniform = (/uniform - /Non-Uniform) //uniform is dis- 
played in Fig. [13] at the three temperatures and is seen to vanish relatively smoothly as 
the transition density nb ~ 0.10 — 0.11 fm~^ is approached. The equivalent pressure differ- 
ence Ap/puniform = (PUniform - PNon-Uniform) /^Uniform and CUtropy difference A///uniform = 

(^Uniform " ^Non-Uniform) /■Suniform, showu in Fig. [14] and Fig. [15] (the latter only at the two 
non-zero temperatures) exhibit discontinuities at the transition density in a way that sug- 
gests the occurrence of a phase transition, as expected. More detailed calculations around 
the transition density are required to determine the exact nature of the transition. Errors 
incurred from the spurious shell energy will have an effect on the details of the results just 
below the transition density and may obscure the nature of the phase transition obtained 
(i.e. first order transition (discontinuity in pressure) or second order (discontinuity in the 
gradient of the pressure)). A more detailed investigation of the phase transition region will 
be presented in a subsequent work. Physically, the transition is expected to be first order, 
as in a liquid-gas transition. However, it can be stated with some confidence that the 3D 
Hartree-Fock model self-consistently predicts a first or second order phase transition from 
non-uniform nuclear matter to uniform nuclear matter. 

D. Dependence on Nuclear Force 

In order to examine the dependence of the energy-deformation surfaces, and ultimately 
the pasta phase diagram, on the properties of uniform nuclear matter (e.g. compressibility, 
symmetry energy), the results of calculations using different nuclear interactions should be 
compared. Fig. [TB] shows plots of the energy-deformation surface at ri]^ = 0.08 fm^^, 7 = 
0° and T = 0.0 MeV as calculated for two different Skyrme parameterizations SkM* and 
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Sly4. The surfaces are qualitatively similar, with the Sly4 force predicting a systematically 
higher free energy-density than SkM* by ~ 10 keV fm~^. Sly4 has a higher symmetry energy 
than SkM* throughout the density and temperature region under investigation, and thus the 
energy cost of adding protons to the matter is higher. A broad minimum in / occurs for both 
parameterizations at A ^ 1000, (3 ~ 0.5 corresponding to a cylindrical hole configuration; 
the free energy-density for that minimum is f — 2.3709 MeV fm~^ for SkM* and / = 2.3816 
MeV fm~^ for Sly4. A more sharply defined minimum occurs for each parameterization 
a,t A — 800, P — 1.0 corresponding to a slab-hke nuclear shape; there, / = 2.3722 MeV 
fm~^ for SkM* and / = 2.3803 MeV fm~^ Sly4. For both Skyrme parameterizations these 
constitute the two lowest energy local minima, although the lower of the two differs between 
the two parameterizations. The difference in energy-density between the minima is ~ 1 keV 
fm~^ in both cases. We can conclude that the two parameterizations offer a qualitatively 
similar spectrum of local minima at this particular density and temperature. A more detailed 
comparison at a variety of densities and temperatures should be carried out and extended to a 
wider range of Skyrme parameterizations (and other compatible energy-density functionals) 
to test the dependence of the pasta phase diagram on properties of the EoS of uniform 
nuclear matter. 



IV. SUMMARY AND CONCLUSIONS 



The self-consistent 3D Skyrme-Hartree-Fock-|-BCS method has been apphed to the study 
of the inhomogeneous 'pasta' phase of bulk nuclear matter through the use of a new code 
TAMAR. A constraint was imposed on both independent quadrupole moments of the neu- 
tron density. Constraining them to a constant value and varying the cell size allows us to 
calculate the variation in free energy- density for a specific nuclear shape as its size (or, cor- 
respondingly, the number of nucleons it comprises) increases. Variation of the quadrupole 
moments then allows a self-consistent determination of the energy-density of different nu- 
clear geometries at a given density and temperature. The variation of all free parameters at 
a given density and temperature plots out an energy-deformation surface / = f{A, /3, 7). In 
this paper the code was applied to the calculation of the properties of inhomogeneous nuclear 
matter in a core-collapse supernova with a fixed value of the proton fraction ?/p = 0.3. 

It has been demonstrated that the present model smoothly spans the transition from iso- 
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lated, roughly spherical nuclei surrounded by a low density gas of nucleons to uniform nuclear 
matter with increasing temperature and baryon number density. The free energy-density of 
non-uniform matter smoothly tends toward that of uniform matter as the transition density 
is reached. Discontinuities in the pressure and entropy are suggestive of a phase transition 
to uniform matter; further calculations are required to determine whether that phase tran- 
sition is first or second order. The five canonical types of pasta are seen to emerge naturally 
from the model as well as a large variety of nuclear shapes intermediate between them. At a 
given density and temperature, many local minima occur corresponding to various nuclear 
geometries, indicating that a variety of pasta shapes coexist; as in neutron star matter [32 1 
the pasta phase in supernova matter is likely to be somewhat disordered. This leads to 
the suggestion that the thermodynamical properties vary smoothly throughout the pasta 
regime rather than a series of phase transitions occurring from one bulk phase made up of 
one specific pasta shape to another: the properties of the different pasta configurations that 
make up the matter are likely to be averaged on a meso- and macroscopic scale. 

Special attention was paid to accounting for all the possible numerical effects that might 
contribute to the results. The choice of original wavefunctions used to start the iteration 
solution to the HF equations and the finite size of the computational cell were shown to 
have no effect on the nuclear geometries obtained. Spurious (numerical) shell effects, arising 
from the artificial discretization of continuum neutron states by the finite computational 
cell have been identified. These effects give rise to oscillatory variation in the energy-density 
of matter with the computational cell size that is distinct from the physical shell effects, 
allowing easy identification when they are the dominant shell energy present. At densities 
just below the density of transition to uniform matter, where the spurious shell effects most 
strongly obscure the physical minima in the energy surfaces, their contribution to the total 
energy-density of matter is in the range ~ 10 — 50 keV fm~^ at A < 400. In order to obtain 
the correct physical energy minima and the transition density with accuracy, the spurious 
shell effects must be taken into account. One possibility is to use analytic expressions for 
the artificial shell energy obtained using the semi-classical WKB approximation 57|]; this 
will be investigated in a subsequent paper. 

The plethora of local minima in the energy-deformation surfaces /(A,/?, 7) arise from 
a combination of physical shell effects from bound nucleons and a number of broad 'ge- 
ometrical' minima occurring over different regions of the energy surface that arise from 
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the competition between the surface energy and Coulomb energies of the different nuclear 
shapes. Physical shell effects that arise from the interaction of the continuum neutrons with 
the nuclear structures, expected to be a significant contribution to the energy spectrum in 



neutron star crustal matter 



32| . play a negligible role in supernova matter as the neutron 



gas external to the nuclear structures is generally an order of magnitude less dense than in 
supernova matter (strongly linked to the fact that the proton fraction is roughly order of 
magnitude higher in supernova matter). Geometrical minima are separated from each other 
by energy barriers in the energy-deformation surfaces of the order 1-5 keV fm~^. At low 
temperatures T < 2.5 MeV, physical shell effects contribute similarly deep local minima. 
The shell effects decrease in magnitude as temperature increases, vanishing aX T ^ 5 MeV 
and above. 

The energy-deformation surface at nb = 0.08 fm~'^, 7 = 0° and T = 0.0 MeV as calculated 
with the SkM* Skyrme parameterization was compared to the same surface calculated by the 
Sly4 parameterization. The surfaces were qualitatively similar, with local minima appearing 
in the same regions of parameter space corresponding to the same nuclear shapes. A more 
detailed examination of the dependence of the pasta phase diagram on the EoS of uniform 
nuclear matter is underway. 

Work is in progress on extending of the present calculation over a larger area of baryon 
number density, proton fraction and temperature in order to construct EoS table for mod- 
eling of core-collapse supernovae and pinpoint more accurately the position in density- 
temperature space and the nature of the phase transition to uniform matter. 
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FIG. 1: (Color) Obtaining a double nuclear shape. The right picture shows a 3D rendering of the 
neutron density profile at nb = 0.08 fm-^, = 0.3, T = 0.0 McV, A = 700, and {f3,j) = (1.0, 
0°). The configuration shown in the left picture is the result of a calculation in a box double the 
size in the z-direction - with A = 1400 - at otherwise the same parameter values. Blue indicates 
the lowest densities and red the highest. 




FIG. 2: Free energy-density / versus nucleon number A aX n\) = 0.11 fm~^, (/3, 7)=(0,0°) and T = 
2.5 MeV. The form of the curve is dominated by a spurious shell energy caused by the finite box 
discretization of the continuum neutron energy spectrum. 
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FIG. 3: Free energy-density / as a function of the number of nucleons A in the unit cell at different 
densities, temperatures and deformations. Top left: (/3,7) = (0.0, 0°), nb = 0.06 fm~^. Top right: 
= (1.0, 0°), nb = 0.06 fm-3. Bottom left: {/3,j) = (0.0, 0°), nb = 0.10 fm-^. Bottom right: 
(/3,^) = (1.0, 0°), nb = 0.10 fm-3. Results at T = 0.0 MeV and T = 2.5 MeV are shown on each 
plot. 
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FIG. 4: Free energy-density / as a function of the number of nucleons A in the unit cell at nb = 0.06 
fm-^ (/3,7) = (0.0, 0°), r = 5 MeV. 
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FIG. 5: (Color) 3D renderings of the neutron density profiles of configurations taken from the 
sequences of increasing A starting at A = 100 and increasing in steps of 100 from left to right and 
from top to bottom (snake fill) up until the last row, where the configurations that deviate from 
the main nuclear geometry are shown (see text for details). Blue indicates the lowest densities and 
red the highest. 
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T = MeV T = 2.5 MeV 




FIG. 6: (Color on-line) Energy-deformation surfaces at nb = 0.08 fm ^, 7 = 60° and T = MeV 
(top left), T = 2.5 MeV (top right), T = 5.0 MeV (bottom left) and T = 7.5 MeV (bottom right) 
as calculated with the Skyrme parameterization SkM*. 




FIG. 7: (Color) 3D renderings of neutron density distributions found at the broad minima on 
the nb = 0.08 fm^^ T = 5.0 MeV; 7=60° energy surface. The left configuration, a cylindrical 
hole, is found at /? = 0.3, A = 800 and the right configuration, a slab-like structure, is found at 
/3 ~ 0.7, A ~ 800. Blue indicates the lowest densities and red the highest. 
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FIG. 8: (Color on-line) Energy-deformation surfaces at rib = 0.06 fm ^, T = 5.0 MeV; 7=0° (left), 
30° (middle) and 60° (right) as calculated with the Skyrme parameterization SkM*. 



nb=0.08 fm"3 nt) = 0.09 fm'^ nb=0.10 fnT^ 




FIG. 9: (Color on-line) Energy-deformation surfaces for T = 2.5 MeV, 7 = 0° and = 0.08 fm ^ 
(left), Ub = 0.09 fm^^ (middle) and nb = 0.10 fm^^ (right). 
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FIG. 10: Neutron density profiles nn(x, y, z) of minimum energy configurations along a line joining 
the centers of adjacent cells in three perpendicular directions x,y and z for the six densities and 
at a temperature of 2.5 MeV. 
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FIG. 11: (Color) 3D renderings of the neutron density profiles of minimum energy configurations 
at T = 2.5 MeV and densities of nb = 0.04 fm~^ (top left), nb = 0.06 fm~^ (top middle), rib = 
0.08 fm^^ (top right), nb = 0.09 fm^'^ (bottom left), rib = 0.10 fm^^ (bottom middle) and nb = 
0.11 fm~^ (bottom right). Blue indicates the lowest densities and red the highest. 



FIG. 12: (Color) 3D renderings of the neutron density profiles of minimum energy configurations 
at nb = 0.10 fm-3 and temperatures of T = MeV (top left), 2.5 MeV (top right), 5 MeV (bottom 
left) and 7.5 MeV (bottom right). Blue indicates the lowest densities and red the highest. 
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FIG. 13: Difference between the free energy-density of uniform matter and non-uniform matter 
A/ divided by the free energy-density of non-uniform matter /uniform as density increases for T = 
MeV (top left), 2.5 MeV (top right) and 5 MeV (bottom). 
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FIG. 14: Difference between the pressure of uniform matter and non-uniform matter Ap divided 
by the pressure of non-uniform matter puniform as density increases for T = MeV (top left), 2.5 
MeV (top right) and 5 MeV (bottom). 
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FIG. 15: Difference between the entropy density of uniform matter and non-uniform matter As 
divided by the entropy density of non-uniforin ms^ttcr suniform ^is density increases for T = 2.5 
MeV (left) and 5 MeV (right). 
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FIG. 16: (Color on-line) Energy-deformation surfaces at n^ = 0.08 fm ^, 7 = 0° and T = 0.0 MeV 
as calculated for two different Skyrme parameterizations SkM* (left) and Sly4 (right). 
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TABLE I: Comparison of nuclear binding energies and root mean square radii of nuclei for three 
different codes using the Sly4 Skyrmc paranictrization. Coulomb and spin-orbit potentials were 
omitted in the calculation. See text for details of the three codes. 



Binding Energy (MeV) Proton (Neutron) RMS Radius (fm) 



Nucleus 


SKYAX 


TDHF3D 


TAMAR 


SKYAX 


TDHF3D 


TAMAR 


160 


126.81 


126.80 


126.80 


2.701 (2.701) 


2.701 (2.701) 


2.701 (2.701) 


40Ca 


399.49 


399.33 


399.35 


3.374 (3.374) 


3.373 (3.373) 


3.373 (3.373) 


56pg 


564.13 


563.96 


563.97 


3.789 (3.875) 


3.789 (3.877) 


3.790 (3.878) 



TABLE II: EoS Quantities vs A for nb = 0.06 fm"^, T = 2.5 MeV, /3,7 = (1.0, 0°). Free energy- 
density / (MeV fm~^), entropy density s (MeV fm~^) and pressure density p (MeV fm"^) are 
given. 



A 


/ 


s 


P 


200 


1.5603 


0.0907 


0.6601 


300 


1.5529 


0.0891 


0.6571 


400 


1.5487 


0.0891 


0.6623 


500 


1.5464 


0.0873 


0.6662 


600 


1.5456 


0.0865 


0.6661 


700 


1.5455 


0.0868 


0.6635 


800 


1.5457 


0.0870 


0.6614 


900 


1.5460 


0.0867 


0.6593 


1000 


1.5464 


0.0862 


0.6573 


1100 


1.5470 


0.0859 


0.6548 


1200 


1.5478 


0.0859 


0.6523 
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TABLE III: EoS at T = MeV. Quantities given are: baryon number density rib (fm"'^), nucleon 
number A and deformation parameters /3, 7 for minimum energy configuration, entropy density 
s (MeV fm~^), free energy-density / (MeV fm~^), pressure p (MeV fm~^), proton, neutron and 
electron chemical potentials /Xp (MeV), (MeV), fie (MeV), and trms divided by half the cell 
width, 6. 





A 


13 


7 


s 


/ 


P 








<5 


0.04 


500 


1.3 


30° 


0.0 


0.901 


0.370 


-32.4 


-1.27 


139.3 


0.71 


0.06 


320 


1.4 


0° 


0.0 


1.593 


0.609 


-35.4 


-1.62 


149.5 


0.88 


0.08 


400 


0.9 


60° 


0.0 


2.364 


0.906 


-37.4 


-1.11 


176.4 


0.88 


0.09 


1000 


0.3 


0° 


0.0 


2.786 


1.019 


-39.6 


-1.80 


182.6 


0.90 


0.10 


1100 


0.0 


0° 


0.0 


3.220 


1.146 


-40.6 


-2.19 


189.2 


0.93 


0.11 


1300 


0.0 


0° 


0.0 


3.662 


1.152 


-42.3 


-3.95 


195.3 


1.00 


0.12 


1200 


0.0 


0° 


0.0 


4.124 


1.447 


-39.7 


-2.59 


201.1 


1.00 



TABLE IV: EoS at T = 2.5 MeV. Quantities as in Table [ml 





A 


P 


7 


s 


/ 


P 








6 


0.04 


400 


1.6 


0° 


0.064 


0.863 


0.393 


-31.9 


-1.89 


139.1 


0.73 


0.06 


340 


1.3 


60° 


0.086 


1.533 


0.654 


-35.5 


-1.90 


159.3 


0.84 


0.08 


800 


0.8 


60° 


0.106 


2.305 


0.947 


-37.9 


-1.77 


175.5 


0.88 


0.09 


1000 


0.2 


0° 


0.115 


2.718 


1.060 


-39.6 


-2.10 


182.5 


0.90 


0.10 


1100 


0.0 


0° 


0.124 


3.145 


1.198 


-41.1 


-2.27 


189.1 


0.93 


0.11 


1300 


0.0 


0° 


0.134 


3.583 


1.299 


-42.7 


-2.84 


195.2 


1.00 


0.12 


1200 


0.0 


0° 


0.137 


4.041 


1.577 


-43.6 


-1.82 


201.0 


1.00 
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TABLE V: EoS at T = 5 MeV. Quantities as in Table [ml 





A 


P 


7 


s 


J 


P 








e 



0.04 


400 


1.6 


0° 


0.261 


0.745 


0.455 


-31.8 


-3.83 


138.6 


0.76 


0.06 


600 


1.1 


60° 


0.356 


1.373 


0.741 


-35.6 


-3.42 


159.0 


0.85 


0.08 


800 


0.4 


0° 


0.437 


2.108 


1.021 


-38.8 


-3.45 


175.1 


0.89 


0.09 


1200 


0.0 


0° 


0.478 


2.501 


1.165 


-40.2 


-3.50 


182.2 


0.92 


0.10 


1200 


0.0 


0° 


0.515 


2.909 


1.327 


-41.6 


-3.41 


188.7 


0.96 


0.11 


1200 


0.0 


0° 


0.545 


3.332 


1.505 


-43.0 


-3.14 


194.9 


1.00 


0.12 


1200 


0.0 


0° 


0.555 


3.784 


1.786 


-43.8 


-1.84 


200.7 


1.00 
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